#!/usr/bin/perl -w
use strict;
use Getopt::Long;
use lib "/net/cpp-group/Leo/bin";
use parse_bl2seq;
use parse_fasta;
#88888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Usage

#88888888888888888888888888888888888888888888888888888888888888888888888888888888
my $usage = <<"USAGE";

USAGE:

    extract_alignment_using_hmm.pl --hmm an.HMM --sequences seq.fa > alignment
		                            [--regex "(\\S+)"]

		This script accepts sequences in the FASTA format and extracts
		the specified domain region
USAGE

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# mandatory


# optional parameters
my $help = undef;
my $hmm = undef;
my $regex = "(\\S+)";
my $seq_file = undef;
GetOptions('help' => \$help,
		   'hmm=s' => \$hmm,
		   'regex=s' => \$regex,
		   'sequences=s' => \$seq_file);

die $usage if ($help);
die $usage unless ($hmm && $seq_file);

my %sequences;
sub save_seq_callback($$$$$$)
{
    my ($acc, $seq, $save_sequences) = @_[0, 1, 5];

	my $sequence = join "", @$seq;
	$sequence =~ tr/BJOUZ//d;
	if ($acc =~ /$regex/)
	{
		$sequences{$1} = $sequence;
	}
	else
	{
		$sequences{$acc} = $sequence;
	}
}
{
	open SEQ, $seq_file or die "Error:\n\tCould not open $seq_file\n$!";
	my $cnt_seq =parse_fasta::parse_sequences(  *SEQ,
												\&save_seq_callback,
												undef);
	close SEQ;
	print STDERR "$cnt_seq\tsequences parsed.\n";

}




open HMMOUTPUT, "/net/cpp-group/tools/hmmer/hmmalign $hmm $seq_file |" or die;

my $cnt_hmm_seq = 0;
while (<HMMOUTPUT>)
{
	last if /#=RF/;
};

chomp;
my ($name, $seq) = split /\s+/;
my %hmm_seq;
$hmm_seq{$name} .=$seq;

while (<HMMOUTPUT>)
{
	chomp;
	($name, $seq) = split /\s+/;
	next unless (defined $seq );
	$hmm_seq{$name} .=$seq;
}

print STDERR scalar keys %hmm_seq, "\tsequences with HMM hits.\n";
die "Error:\n\tHMM Reference sequence absent\n" unless (exists $hmm_seq{'#=RF'});
$hmm_seq{'#=RF'} =~ /^(\.*)(x[x.]*?)(\.*)$/
 or die "Error:\n\tHMM Reference sequence in the wrong format\n";
my $start	  = length($1);
my $align_len = length($2);
#print STDERR "$start\n$align_len\n";
my %seq_gaps;
for my $acc(keys %hmm_seq)
{
	next if ($acc eq '#=RF');
	my $seq = $hmm_seq{$acc};
	my @gap_str;
	while ($seq =~ /^([^.\-]*)([.\-]+)(.*)$/)
	{
		push(@gap_str, length($1).":".length($2));
		$seq = $1.$3;
	}
	my $gap_str = join ("|", @gap_str);
	#print STDERR ">$acc\n$hmm_seq{$acc}\n";
	#print STDERR expand_prot_with_alignment_str($seq, '-', $gap_str), "\n";
	#exit;
	$seq_gaps{$acc} = $gap_str;
}

my @acc = keys %seq_gaps;
for (my $i = 1; $i < @acc; ++$i)
{
	for (my $j = 0; $j < $i; ++$j)
	{
		print join ("\t",
						$acc[$i],
						$acc[$j],
						$start,
						$align_len,
						$align_len,
						1.0,
						0,
						length($hmm_seq{$acc[$i]}),
						length($hmm_seq{$acc[$j]}),
						$align_len,
						0,
						$seq_gaps{$acc[$i]},
						$seq_gaps{$acc[$j]}), "\n";
	}
}

#seq1	seq2	start	len	positions	coverage	e_value	len1	len2	residues	bitscore	gaps1	gaps2

